home *** CD-ROM | disk | FTP | other *** search
/ Info-Mac 4 / Info_Mac IV CD-ROM (Pacific HiTech Inc.)(August 1994).iso / Science / RLaB / examples / house.r < prev    next >
Text File  |  1994-04-25  |  2KB  |  108 lines
  1. //
  2. // From MATRIX Computations, G.H. Golub, C.F. Van Loan
  3. //
  4.  
  5. //
  6. // Create a LIST to hold the Householder related functions.
  7. // These functions have only been used with REAL inputs.
  8. //
  9.  
  10. house = <<>>;
  11.  
  12. //
  13. // house.v(): Given an N-vector X, generate an n-vector V
  14. // with V[1] = 1, such that (eye(n,n) - 2*(V*V')/(V'*V))*X
  15. // is zero in all but the 1st component.
  16. //
  17.  
  18.   house.v = function(x)
  19.   {
  20.     local(b, n, u, v);
  21.   
  22.     v = x;
  23.     n = max( size(x) );
  24.     u = norm(x, "2");
  25.     if(u != 0)
  26.       {
  27.     b = x[1] + sign(x[1])*u;
  28.     if(n > 1) {
  29.       v[2:n] = v[2:n]/b;
  30.     }
  31.       }
  32.     v[1] = 1;
  33.     return v;
  34.   };
  35.  
  36. //
  37. // house.row(): Given an MxN matrix A and a non-zero M-vector V
  38. // with V[1] = 1, the following algorithm overwrites A with 
  39. // P*A, where P = eye(m,m) - 2*(V*V')/(V'*V)
  40. //
  41.  
  42.   house.row = function(A, v)
  43.   {
  44.     local(a, b, w);
  45.     
  46.     a = A;
  47.     b = -2/(v'*v);
  48.     w = b*A'*v;
  49.     a = A + v*w';
  50.     return a;
  51.   };
  52.  
  53. //
  54. // house.col(): Given an MxN matrix A, and an N-vector V, 
  55. // with V[1] = 1, the following algorithm overwrites A with A*P
  56. // where P = eye(N,N) - 2*(V*V')/(V'*V)
  57. //
  58.  
  59.   house.col = function(A,v)
  60.   {
  61.     local(a, v, b, w);
  62.  
  63.     b = -2/(v'*v);
  64.     w = b*A*v;
  65.     a = A + w*v';
  66.  
  67.     return a;
  68.   };
  69.  
  70. //
  71. // Given A, with M >= N, the following function finds Householder
  72. // matrices H1,...Hn, such that if Q = H1*...Hn, then Q'*A = R is
  73. // upper triangular.
  74.  
  75. // House.qr returns a MxN matrix, with the upper triangular part 
  76. // containing [R]
  77.  
  78. house.qr = function(A)
  79.   {
  80.     local(a,j,n,m,v);
  81.  
  82.     a = A;
  83.     m = size(A)[1]; n = size(A)[2];
  84.     v = zeros(m,1);
  85.  
  86.     for(j in 1:n)
  87.       {
  88.         v[j:m] = house.v( a[j:m;j] );
  89.         a[j:m;j:n] = house.row( a[j:m;j:n], v[j:m] );
  90.     if(j < m) {
  91.     a[ (j+1):m;j ] = v[(j+1):m];
  92.         }
  93.       }
  94.     return a;
  95.   };
  96.  
  97. //
  98. //  Generate P matrix
  99. //
  100.  
  101. P = function(V)
  102. {
  103.   local( m );
  104.  
  105.   m = max( size(V) );
  106.   return [ eye(m,m) - 2*(V*V')./(V'*V) ];
  107. };
  108.